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■ We present the results of a large library of cosmological A^-body simulations, using 
psj I power-law initial spectra. The nonlinear evolution of the matter power spectra is com- 
^3 ' pared with the predictions of existing analytic scaling formulae based on the work 

of Hamilton et al. The scaling approach has assumed that highly nonlinear struc- 
tures obey 'stable clustering' and are frozen in proper coordinates. Our results show 
that, when transformed under the self-similarity scaling, the scale-free spectra de- 
O ' fine a nonlinear locus that is clearly shallower than would be required under stable 

clustering. Furthermore, the small-scale nonlinear power increases as both the power 
22 ' spectrum index n and the density parameter O decrease, and this evolution is not 

well accounted for by the previous scaling formulae. This breakdown of stable clus- 
tering can be understood as resulting from the modification of dark-matter haloes by 
continuing mergers. These effects are naturally included in the analytic 'halo model' 
. , for nonlinear structure; we use this approach to fit both our scale-free results and 

■ also our previous CDM data. This method is more accurate than the commonly- 
used Peacock-Dodds formula and should be applicable to more general power spec- 
tra. Code to evaluate nonlinear power spectra using this method is available from 
http : //asl . chem.nottingham. ac .uk/~res/sof tware .html. Following publication, 
we will make the power-law simulation data publically available through the Virgo 
website http: //www. mpa-garching.mpg.de/Virgo/. 
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ABSTRACT 



o 



1 INTRODUCTION 

In the current cosmological paradigm, structures grow 
through the gravitational instability of coUisionless dark 
matter fluctuations. This occurs in a hierarchical way, with 
small-scale perturbations collapsing first and large-scale per- 
turbations later. One of the most direct manifestations of 
this nonlinear process is the evolution of the power spectrum 
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of the mass, P{k), where k is the wavenumber of a given 
Fourier mode. Understanding this evolution of the power 
spectrum is one of the key problems in structure formation, 
being directly related to the abundance and clustering of 
galaxy systems as a function of mass and redshift. If the 
processes that contribute to the evolution can be captured 
in an accurate analytic model, this opens the way to us- 
ing observations of the nonlinear mass distribution (from 
large-scale galaxy clustering or weak gravitational lensing) 
in order to recover the primordial spectrum of fluctuations. 
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One of the most influential attempts at such an analytic 
description of clustering evolution was the 'scaling ansatz' 
of Hamilton et al. (1991; HKLM), which is described in Sec- 
tion 2. This scaling procedure was generalized to models 
with Q ^ 1 and given a more accurate A'^-body calibration 
by Peacock & Dodds (1996; PD96). HKLM assumed that a 
nonlinear collapsed object would decouple from the global 
expansion of the Universe to form an isolated system in virial 
equilibrium - the 'stable clustering' hypothesis of Davis & 
Peebles (1977). This assumption has boon widely adopted, 
and yet it appears somewhat inconsistent with hierarchical 
models - in which objects are continuously accreting mass 
and growing through mergers. Indeed, the validity of stable 
clustering has been increasingly questioned in recent years 
(e.g. Yano & Gouda 2000; Caldwell ot al. 2001). One of our 
aims in this paper is thus to establish whether stable clus- 
tering is relevant for understanding the small-scale evolution 
of the power spectrum. 

We therefore explore the gravitational instability of 
dark matter fluctuations through a scries of large A'^-body 
simulations of clustering from power-law initial conditions, 
with 



P{k) oc k" 



(1) 



We consider both n = 1 models, in which the evolution 
can obey a similarity solution, and also low-density models 
with and without a cosmological constant. We demonstrate 
that the resolution of the simulations is sufliciont to mea- 
sure the power well into the regime at which the HKLM 
procedure predicts a well-defined slope for the power spec- 
trum determined by stable clustering. In practice, we find 
that the power spectra are generally shallower than would 
be required for clustering to bo stable on small scales. Fur- 
thermore, as both n and f2 decrease, the amplitude of the 
small-scale spectrum increases in a maimer that is not well 
described by any of the previous fitting formulae. In light of 
these results, a new method for predicting nonlinear spectra 
is proposed. This method is based on the 'halo model' (e.g. 
Seljak 2000; Peacock & Smith 2000), which does not assume 
stable clustering. This allows us to fit our data and also the 
cold dark matter (CDM) data of Jenkins et al. (1998; J98) 
with a high degree of accuracy. 

The paper is structured as follows. In Section 2 we pro- 
vide a brief overview of the theoretical understanding of 
nonlinear evolution. In particular, a description of the stable 
clustering hypothesis, the nonlinear HKLM scaling relations 
and the halo model are given, as these ideas are central to 
this paper. We also discuss the scale-free models and their 
self-similarity properties. In Section 3 we describe the nu- 
merical simulations and we provide a visual comparison of 
the growth of structure in the different scale-free models. 
In Section 4 we describe an improved method for measur- 
ing power spectra and in Section 5 we present the power 
spectra data and contrast them with the current nonlinear 
fitting formulae. In Section 6 we describe a new approach 
to fitting power spectra and its generalization to CDM, and 
then compare our new globally optimized formula with the 
results from Section 5 and also the CDM data. Finally, in 
Section 7 we draw our conclusions and discuss our findings 
in a wider context. 



2 DESCRIPTION OF NONLINEAR 
EVOLUTION 

2.1 Prom linear theory to stable clustering 

The mass density field, at comoving position x and time t, 

is defined as 



p(x,t)=p(t)[l-F5(x,i)] , 



(2) 



where 5 is the density fluctuation about the homogeneous 
background p. The 2-point auto-correlation function of the 

density field is 

C(r) = (<5(x)<5(x + r)) , (3) 

which in three dimensions is related to the dimensionless 
power spectrum A^(/c) through the integral relation 

sinfcr dk 



ar) = / A^(fc) 



kr k ' 



(4) 



where we have assumed that the field is isotropic and ho- 
mogeneous, is the contribution to the firactional density 
variance per unit In k. In the convention of Peebles (1980), 

this is 



A\k) 



V 



dink (27r)- 



■4nk^P{k) 



(5) 



V being a normalization volume. 

When (5(x, t) <C 1 the temporal evolution of the fluctu- 
ation is separable and the field scales as 

S{^,t) = ^^S{^,to) , (6) 

where D{t) is a growth factor whose exact form can be deter- 
mined from linear theory. As (5(x,t) 1, increasingly higher 
orders of perturbation theory are required (see Bernardeau 
et al. 2001 for a thorough review). Eventually, perturbation 
theory fails and numerical methods must be applied. Even 
so, it was proposed (Peebles 1974a; Davis & Peebles 1977; 
Peebles 1980) that clustering in the very nonlinear regime 
might be understood by assuming that regions of high den- 
sity contrast undergo virialization and subsequently main- 
tain a fixed proper density. The correlation function for a 
population of such systems would then simply evolve accord- 
ing to ^(r, t) oc 1/p (X a^, where r is a proper distance. This 
evolution was termed 'stable clustering'. Peebles went on to 
show that if the initial power spectrum was a pure power- 
law in k with spectral index n, P{k) cx k", and 'd — 1, 
then under the stable clustering hypothesis, the slope of the 
nonlinear correlation function would be directly related to 
the spectral index through the relation 



^{r,t) oc r 



7 = 



3(3 -I- n) 



(7) 



5 + n 

Hence, if stable clustering applies, then the nonlinear density 
field retains some memory of its initial configuration, and in 
principle can be used to measure the primordial spectrum 
of fluctuations. 

2.2 The HKLM scaling relations 

HKLM developed a method for interpolating between linear 
theory on large scales and the nonlinear predictions of the 
stable clustering hypothesis on small scales. They showed 
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that the nonhnear volume averaged two-point correlation 
function, 



^(^) = j3 / y ^^y^ ^y > 



(8) 



measured from the scale- free simulations of Efstathiou et al. 
(1988), could be parameterized by a simple function of the 
linear correlation function, provided that nonlinear evolu- 
tion were to induce a change of scale. 

The transformation of scales follows from an intuitive 
continuity argument, based upon the 'spherical top-hat' 
model. Let the mass enclosed within a spherical overden- 
sity in the initial stages of evolution be mo (< i) and its 
mass at some later time be m (< x). As each shell evolves, 
it will reach a maximum expansion point, turn around and 
collapse. If there is no shell crossing, then mass is conserved 
and 

4 q 4 o 

mo (< £) = -7rp (< £) t = -ivp (< x) x= m (< a;)) . (9) 

o o 

The argument now identifies 1 -|- |^ as the factor by which 
the density is enhanced relative to the mean (Peebles 1980). 
Provided ^ 1, this implies the scaling 



x"" [l+CNL(a;,t)] 



(10) 



where x represents a nonlinear scale and £ a Lagrangian 
scale. 

Finally, after this rescaling, the nonlinear correlations 
are taken to be a universal function of the linear ones: 



^Mx,t) = f [^r.{£,t)] 



(11) 



HKLM then assumed that the functional form of f{y) could 
be determined analytically in two regimes: in the linear 
regime, where <C 1, ,f{y) = y ; when ^ 1, galaxy 
groups would exhibit 'stable clustering', for which A^l a" 
and since Al oc a^, this implied that f{y) oc t/^^^. The 
interpolation between these two regimes, where j/ ~ 1, 
was determined empirically by HKLM, by comparison with 
numerical simulation. However, Padmanabhan (1996) pro- 
posed that the quasilinear regime could also be understood 
analytically. He considered the point at which a spheri- 
cal perturbation would reach its maximum radius, which 
is Xmax = 1/Sl oc l/^L, according to the spherical model. 
Padmanabhan thus conjectured that 



(12) 



f . m mo I 

oc p(< a;^a^) oc _ oc -3— oc 

(in effect rediscovering the argument of Gott & Roes 1975). 
Although useful heuristically in explaining why the quasilin- 
ear regime of /nl should be steeper than either the linear or 
nonlinear regime, it is not clear that this expression matches 
the observed quasilinear slope very well (Padmanabhan et 
al. 1996; Jainl997). Wc investigate this further in Section 5. 

HKLM's nonlinear scaling argument was further devel- 
oped by Peacock & Dodds (1994; PD94), who proposed that 
the scaling ansatz could be used for predicting power spectra 
by simply replacing ^ — > A^ and letting the linear and non- 
linear scales represent linear and nonlinear wavcnumbers: 
£ = k^^ and x — k^l- This suggested the formalism 



.(^nl) =/nl [Ai(fcO] 



The accuracy of the HKLM and PD94 scaling formulae 
was tested by Jain, Mo & White (1995; JMW95). They per- 
formed a series of simulations with 100^ particles as opposed 
to the previous 32^, and discovered that the nonlinear lo- 
cus described by the data exhibited a strong n-dcpendencc. 
The HKLM and PD94 functions underestimated the mea- 
sured correlation functions and power spectra, the fits being 
worse for more negative n. JMW95 then showed that this 
n-dependence could be removed by a simple scaling of the 
variables in the log^NL(a;,t) — log^L(^, i) plane. In order for 
the model to be applied to curved spectra, such as the CDM 
model, an effective spectral index riefi was required. JMW95 
proposed that the appropriate n should be given by 



ries 



dlnP{k) 



dink 



(14) 



where Rc is the scale on which the variance of the density 
field is unity. This showed the right response with scale, and 
described their data to a precision of 15 — 20%, which was 
adequate given the scatter within the simulations. 

Further refinements were again made by Peacock & 
Dodds (1996; PD96), who used a large ensemble of 80'^ par- 
ticle simulations to investigate the n-dependence and the 
response of the clustering to low density universes: f2 < 1 
and f2 -h A = 1, where Q and A are the densities associated 
with matter and the cosmological constant, relative to the 
Einstein-de Sitter universe. PD96 concluded that nonlinear 
effects tend to increase the power on small scales for spectra 
with more negative spectral indices and for lower densities. 
PD96 also produced a fitting formula which modelled their 
data, and also CDM-like spectra through defining an effec- 
tive spectral index that changed with each wavenumber 



, , d In P 
ncfi(fcL) = -r. — r(fc = Kl/2) 
am K 



(15) 



fcNL = [1 + Anl(A:nl)] k 



(13) 



Subsequently, high resolution numerical simulations of 
CDM-like universes have shown that the PD96 formulae 
match the observed nonlinear power spectra closely (Mo, 
Jing & Borner 1997; J98; Smith et al. 1998), but with some 
significant deviations. Jain & Bcrtschinger (1998) found a 
larger discrepancy in their 256^ P^M simulation of cluster- 
ing from an n = — 2 power spectrum, with both the for- 
mula of JMW95 and PD96 underestimating the quasi-linear 
power. They also claimed that their results for highly non- 
linear clustering were in accordance with stable clustering, 
although finite volume effects have drawn their results into 
question (Ma & Fry 2000a; Scoccimarro et al. 2001). We 
discuss this issue in further detail in Section 3.3. Recent 
attempts to constrain cosmological parameters from weak 
gravitational lensing studies, that require as input the non- 
linear matter power spectrum, have also uncovered deficien- 
cies in the PD96 formula, with the poorest performance for 
the a = 1 tCDM model (Van Waerbeke et al. 2001). 

2.3 A dark matter halo approach 

More recently an entirely different analytical model for non- 
linear gravitational clustering has emerged: the 'halo model'. 
In this model, the density field is decomposed into a distri- 
bution of clumps of matter with some density profile. This 
basic idea goes back to Neyman & Scott (1952), and recurs 
in more modern form in (Scherrer & Bcrtschinger 1991). 
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Following the realization that galaxy bias was strongly in- 
fluenced by the number of galaxies in a halo (Jing, Mo & 
Borner 1998; Benson et al. 2000), a number of authors (Sel- 
jak 2000; Peacock & Smith 2000; Ma & Pry 2000a; Scoc- 
cimarro ct al. 2001) resurrected the Neyman-Scott model 
with a modern mass function for dark haloes (Press & 
Schechter 1974; Sheth & Tormen 1999; Jenkins et al. 2001), 
plus realistic density profiles (Navarro, Frenk & White 1996; 
Navarro, Prenk & White 1997; Moore et al. 1999), and a 
mass-dopcndcnt galaxy 'occupation number'. The inclusion 
of bias is an attractive aspect of the halo model, but we will 
not be concerned with this here. 

In the halo model, the large-scale clustering of the mass 
arises through the correlations between different haloes. Pre- 
scriptions for this clustering wore given by Mo & White 
(1996); Mo, Jing & White (1997); Sheth & Lemson (1999); 
Sheth & Tormen (1999); Sheth, Mo & Tormen (2000), and 
a recent example of their effectiveness is shown clearly in 
Colberg et al. (2000). On small scales, the correlations are 
derived purely from the convolution of the density profile of 
the halo with itself (Peebles 1974b; McClelland & Silk 1977; 
Sheth & Jain 1997). This model thus makes strong predic- 
tions about the clustering on small scales. Unless the den- 
sity profile and mass function obey a specific relationship, 
the merger-driven evolution of the mass function means that 
stable clustering approximation does not hold true (Yano & 
Gouda 2000; Ma & Pry 2000b). For a more detailed review 
of the halo model and its applications we refer the reader to 
Sheth & Cooray (2002). 



2.4 Scale-free models 

An elegant way to study nonlinear evolution is to simulate 
'scale free' universes that have no inbuilt characteristic phys- 
ical length scales. We follow Efstathiou et al. (1988) and 
require two conditions to be satisfied: 

(1) The initial power spectrum of fluctuations is a power 
law: 



P{k) = Ak" 



1 < n < -3. 



(16) 



(2) The evolution of the scale factor for the cosmological 
model power law in time: 



a{t) oc t° 



(17) 



The most interesting cosmological model that satisfies these 
constraints is the Einstein-de Sitter model: a = 2/3, n = 1 
and A = 0, so that the linear-theory growth of the power 
spectrum is P{k) oc a^. 

In this case, the only natural way to define a charac- 
teristic length is through the scale at which the fluctuations 
become nonlinear. The variance of the linear density fleld, 
smoothed on some comoving length scale x, is 



a^(x,a) = j 



Al{k,a) \Wikx)f 



k 



(18) 



where W is the filter function. If we assume A^(o, fe) oc 
o^fe^"*"", and that the filter causes a cut-off at some high 
spatial frequency fcc ~ 1/a;, we find 



a\kc,a) oc / a'^k'^+^dk oc a'x:^^+''^ . 
Jo 



We now define a nonlinear wavenumber, A;nl such that 
iT^(feNL,a) = 1, so that 



A;NL(a)oca-^/(^+"). 



(20) 



Under this transformation, it is plausible that the statistics 
of gravitational clustering will be expressible as a similarity 
solution: 



P(fc,a) = P(fc/fcNL) 



(21) 



(19) 



(Davis & Peebles 1977; Peebles 1980; Efstathiou et al. 1988; 
Jain & Berchinger 1998). No formal proof of the similarity 
solution exists, and this conjecture is something that must 
be tested empirically via simulation. We refer the reader 
to the work of Colombi, Bouchet & Hernquist (1996) for 
further discussion of the range of spectral indices for which 
self-similarity should be valid. 

In practice, we present good evidence in this paper that 
the power spectrum does scale in this way for > n > —2. 
Spectra outside this range are harder to simulate and so 
not yet tested. We may however anticipate that only cer- 
tain initial spectra will evolve in a self-similar fashion. For 
n > 1, the amplitude of gravitational potential fluctuations 
diverges on small scales, so one might question the idea 
of a hierarchy that grows via the merger and disruption 
of small systems. However, this argument is not definitive, 
since the similarity solutions generally depart from P oc fc" 
for k > Ainl. We seek a function which is of this power- 
law form for k < /cnl and some unknown form at larger fc, 
and which evolves in a self similar fashion. In practice, this 
function is found by starting with exact power-law initial 
conditions, and hoping that the simulation will relax into 
the desired self-similar form as it evolves. The existence of 
a self-similar solution with n > 1 on large scales therefore 
remains an open question. On large scales, the peculiar ve- 
locity field diverges if n < — 1, so more negative indices may 
seem problematic. This does not seem to be a problem in 
practice, probably for the reasons discussed by Bernardeau 
et al. (2001): the divergent modes of very long wavelength re- 
ally just cause a translation, and Galilean invariance means 
that the statistics of smaller-scale clustering are unaffected. 
Certainly, well-defined results can be obtained from pertur- 
bation theory for n more negative than —1, so the only clear 
limit is n < — 3, for which the whole idea of asymptotic ho- 
mogeneity breaks down. 

If we can find initial spectra for which self-similarity ap- 
plies, this is an extremely useful means of assessing the reli- 
ability of iV-body results. Also, over limited ranges of mass, 
the scale-free models correspond directly to more physically 
motivated models such as CDM, whose spectral index is a 
slow function of scale. As we shall show, an analytic de- 
scription of nonlinear evolution in the scale-free case leads 
quite directly to a method that can also give an accurate 
description of nonlinear evolution in CDM models. 



3 THE NUMERICAL SIMULATIONS 

We have produced a large library of A'^-body simulations 
with N — 256'^ particles. We considered Einstein-de Sitter 
(O = 1) models, and also low-density open and fiat A ge- 
ometries. The spectral indices that have been simulated are 
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Table 1. Parameters of the 256"^ particle, scale-free simulations. The rl simulations represent glass initial conditions and r2 simulations 
are grid starts. 



simulation 


e/L 


A2(fcb,a = 1) 


^initial 


«flnal 


timesteps 


energy error 


output values of a 


n = 


-2 rl 


0.00025 


0.133 


0.025 


0.62 


831 


0.04 % 


0.025, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.62 


n = 


-2 r2 


0.00025 


0.133 


0.025 


0.55 


904 


0.04 % 


0.025, 0.1, 0.2, 0.3, 0.4, 0.5, 0.55 


n = 


-1.5 rl 


0.00023 


0.046 


0.010 


0.96 


991 


0.16 % 


0.01, 0.25, 0.315, 0.4, 0.5, 0.63, 0.794, 0.96 


n = 


-1.5 r2 


0.00023 


0.046 


0.010 


1.00 


915 


0.16 % 


0.01, 0.25, 0.315, 0.4, 0.5, 0.63, 0.794, 1.0 


n = 


-1 rl 


0.00023 


0.017 


0.010 


0.83 


991 


0.31 % 


0.01, 0.25, 0.315, 0.4, 0.5, 0.63, 0.794, 0.83 


n = 


-1 r2 


0.00023 


0.017 


0.010 


1.00 


815 


0.31 % 


0.01, 0.25, 0.315, 0.4, 0.5, 0.63, 0.794, 1.0 


n = 


rl 


0.00025 


0.003 


0.025 


0.66 


1443 


0.50 % 


0.025, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.66 


n = 


r2 


0.00025 


0.003 


0.025 


0.50 


1239 


0.50 % 


0.025,0.1,0.2,0.3,0.4,0.5 



n = —2, —1.5, —1 and and two realizations of each spec- 
tral index were carried out. The simulations were executed 
on either 128 or 64 processors of the Edinburgh Cray T3E 
supercomputer, using the parallelized P^M 'shmem' version 
of HYDRA (Macfarland ct al. 1998; Couchman, Thomas & 
Pearce 1995; Pearce & Couchman 1997), in purely collision- 
less dark matter mode. 

The large-scale force calculation in HYDRA used a 512^ 
Fourier mesh, supplemented by direct summation of close 
pairs to achieve the desire total interparticle force. As usual, 
this is softened on small scales in order to suppress two-body 
encounters. In HYDRA, the transition from pure Newtonian 
to constant force is achieved using a 'spline- kernel softening': 
with this method, the interparticle forces become precisely 
Newtonian after 2.34 times the softening length. In all cases, 
we adopted a comoving softening length that is simply a 
fraction / of the interparticle spacing 

£ = , (22) 

where L is the side of the simulation box. We used / ~ 0.064, 

which is slightly smaller than the late-time value used by 
Efstathiou et al. (1988) and the small-box calculations of 
J98 who used / ~ 0.1. However, it is slightly larger than 
the values used by Jain & Bertschinger (1998) who used an 
effective value of / = 0.05, and also the value chosen by 
J98 for their big-box simulations, / ~ 0.038. We ran a few 
test simulations in which / was varied, and we believe that 
the results quoted here are not sensitive to the exact value 
adopted. 

For the initial particle load, a combination of 'quiet' 
starts and 'glass' configurations was used. The quiet starts 
were produced by simply placing particles onto a uniform 
grid with spacing L/N^^^ . This method gives no contribu- 
tion to the power spectrum from particle placement except 
on scales of the order half a mesh spacing (See Section 4). 
However, grid initial conditions may lead to non-physical 
features on very small scales at late times. An example of 
this occurs in the Warm Dark Matter simulations of Bode, 
Ostriker & Turok (2001), where the population of 'secondary 
objects' which they find to form by fragmentation of sheets 
and filaments may actually be a rmmerical artefact induced 
by the grid. An alternative approach is the glass-like dis- 
tribution that is obtained when a random distribution of 
particles is evolved with the signs of the N-body accelera- 
tions reversed (White 1993; Baugh, Gaztanaga & Efstathiou 
1995). The resultant particle distribution displays no regu- 
lar pattern, but is sub-random. By construction, the glass 



initial conditions are non-evolving in the absence of pertur- 
bations. The glass load was generated once, but can be used 
in many different simulations by adding in the appropri- 
ate displacement field. This was generated from the initial 
density field using the approximation of Zel'dovich (1970). 
The Fourier modes of the density field were a Gaussian re- 
alization, with random phases and amplitudes chosen from 
a Rayleigh distribution. 

For both the grid and glass methods, particle discrete- 
ness on the smallest scales leads to a spectrum that is com- 
parable to that of the shot-noise distribution on that scale. 
Numerical evolution should proceed until the scales of in- 
terest are well above this noise. For most spectra, memory 
of the initial small-scale discreteness is only truly lost after 
expansion by roughly a factor of 10 (see Section 4.2). 

3.1 Self-similar simulations 

The normalization of the scale-free power spectra is most 
simply specified in terms of the power on the box scale at 
the epoch when the expansion factor a is unity, 

/£) ^^^^ 

where = 2'n/L. The benefit of normalizing the spectrum 

in this way is that the box-scale power is directly related to 
the error induced through omitting modes with wavelength 
above L, and so the effects can be monitored (see Section 
3.3). 

Table 1 displays all relevant simulation parameters 
for the scale-free runs. A large degree of nonlinearity was 
achieved for all of the simulations and the n — —1 and —1.5 
calculations were completed to the specified level of normal- 
ization. The n = calculations were halted after the cube 
had expanded by roughly a factor of 25, due to the intense 
demands on the cpu time from performing the PP part of 
the calculation. Also, the n = —2 calculations were halted 
after a similar factor of growth; this was due to the problems 
of finite volume effects, which we discuss in detail in Section 
3.3. 

Figs 1 and 2 provide a visual account of the growth of 
structure in the four models. We show three epochs from the 
four different models: the initial conditions; an intermediate 
epoch; the final output epoch. The n — —2 simulations dis- 
play a number of large-scale fiuctuations which collapse to 
form large filaments and groups, whereas the n = simula- 
tions are characterized by a large number of tightly bound 
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Figure 1. Slices siiowing the growtli of structure in the glass n = —2 simulation (left column) and 'grid-start' n = —1.5 simulation (right 
column). All of the slices are of thickness L/10. Prom the n = —2 simulation we show expansion factors a = 0.2, 0.45 and 0.55, and from 
the n = —1.5 simulation we show epochs a = 0.25,0.63 and 1.0. The normalization of the final states in the n = —1 and n = —1.5 runs 
were A^{2n/L,a = 1.0) = 0.133 and 0.046, respectively. 
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Figure 2. Same as Fig. 1, but this time sliowing the comoving projection of particles in the glass n = —1 simulation (left column) and 
glass n = simulation (right column). Prom the n = —1.0 simulation we show epochs a = 0.25,0.63 and 0.83, and from the n = 
simulation we show expansion factors a = 0.1,0.3 and 0.5. The normalization of the final states in the n = —1 and n = runs were 
A2(27r/L, a = 1.0) = 0.017 and 0.003, respectively. 
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Table 2. Parameters of the 256^ particle, power-law A and open simulations. Epochs include a = 0.025, 0.05, 0.1, 0.2, 0.25, 0.3, 0.4, 0.5, 
0.6, 0.7, 0.8, 0.9, 1.0 



simulation 


e/L 


A2(fcb) 


n 


A 


^initial 


Iflnal 


timcsteps 


energy error 


n = -2 


0.00025 


0.0479 


0.26 


0.74 


0.025 


1.0 


1065 


0.05 % 


n = -2 


0.00025 


0.0479 


0.2 


0.0 


0.025 


1.0 


965 


0.09 % 


n = -1.5 


0.00025 


0.0240 


0.26 


0.74 


0.025 


1.0 


971 


0.13 % 


n = -1.5 


0.00025 


0.0240 


0.2 


0.0 


0.025 


1.0 


965 


0.13 % 


n = -1 


0.00025 


0.0101 


0.26 


0.74 


0.025 


1.0 


1342 


0.28 % 


n = -1 


0.00025 


0.0101 


0.2 


0.0 


0.025 


1.0 


965 


0.87 % 


71 = 


0.00025 


0.0003 


0.26 


0.74 


0.025 


1.0 


1020 


0.86 % 


n = 


0.00025 


0.0003 


0.2 


0.0 


0.025 


1.0 


967 


1.86 % 



objects and a paucity of large-scale filamentaxy features, in 
accordance with the results of Efstathiou et al. (1988). Fig. 

1 also compares the glass start to the grid starts. In the glass 
start no features other than the prescribed fluctuations are 
observed, whereas the grid start shows faint lattice patterns 
which are still observable in the voids at the final epoch. 



3.2 Power-law open and flat simulations 

At late times the amplitude of the nonlinear power spectrum 
is very sensitive to the density of the universe, and strongly 
modulates the amplitude of the nonlinear clustering signal. 
This effect is important to quantify if one wishes to con- 
struct general models for evolving nonlinear power spectra. 
We investigated this density dependence by performing a 
further series of high resolution, 256"' particle, simulations 
for open universes where ft = 0.2 at the final epoch and 
for flat universes where Q = 0.26 and A = 0.74 at the final 
epoch. The values for the density parameter were selected 
so that each full integration would span a large dynamic 
range in fl. The amplitude of the final box-scale mode was 
set slightly lower than in the f2 = 1 simulations, because 
of the greater small-scale nonlinearities that are generated 
in low-density models. For all of these simulations wc have 
used the glass initial particle load. Table 2 displays all of the 
relevant simulation parameters. 



3.3 The challenge of n ^ -3 

On small scales, the slope of the CDM power spectrum ap- 
proaches n ~ —3, so it is important to understand how such 
spectra evolve in the nonlinear regime. However, highly neg- 
ative spectral indices have proven difficult to simulate (Efs- 
tathiou et al. 1988; Jain, Mo & White 1995; PD96; Jain & 
Bertschinger 1998), and this can be attributed to two main 
effects. 

First, the number of particles must be high enough to 
simulate virialized clusters convincingly. Second, the finite 
size of the simulation volume means that the longest wave- 
length fluctuations that are present arc Xh = L ; kh = 2tt/L. 
The absence of modes beyond the box scale induces an error 
in the nonlinear spectrum, since nonlinearity couples Fourier 
modes together and power leaks from large to small scales; 
the importance of this effect increases for increasingly neg- 
ative spectral indices and dominates as n approaches —3. 

The error in the power spectrum due to these missing 



modes can be estimated from the linear power spectrum. We 
can quantify the missing variance as follows: 

- 1 Ai(fc') 



Al(fc) 



dk \ ^ 1 

T ^ ^ 4^(^2-hm2-hn2)3/2 



(24) 



where the sum is over all integer triples {£, m, n) except 
(0,0,0) and the wavenumber k' = kh{f + + n^Y''^ . 
Strictly speaking, both terms on the rhs are divergent for 
power-law spectra with n > —3. Nonetheless, if one imposes 
a sufficiently smooth cut-off at fecut in the power spectrum, 
then the difference is well defined in the limit of fecut oo. 

We have estimated (T^iiss numerically in this way for 
scale-free power spectra as a function of n. To about 1% 
accuracy the result is given by: 

2 _ Ai(fcb 



3-l-n 



■F(3-hn), 



(25) 



where F(y) = 1-0. 31y+0.015y^ -1-0. 001332/^ and this expres- 
sion is valid for —3 < n < 1. One can check the numerical 
result, not only by confirming it is insensitive to the precise 
value of fccut, but also for the special case n = where it is 
easy to see from geometric considerations that the value of 
F(3) is 3/47r. In the limit of n ^ —3 the missing variance is 
well approximated by the quantity cj\^^ defined as: 

(26) 

So as to ensure that the missing variance does not become 
significant for our simulations, we have chosen to adopt the 
criterion 



o-L < 0.04 , 



(27) 



for which the large-scale missing modes are safely linear. 

It is for these reasons that the relatively low resolution 
(compared to modern standards) 32^ particle, n = — 2 simu- 
lation of Efstathiou et al. (1988) could only reproduce the ex- 
act similarity solution for the power spectrum over a narrow 
range of expansion. Also, for the more recent high-resolution 
256^ particle, n = — 2 simulation of Jain & Bertschinger 
(1998), the box-scale power for their last three outputs vi- 
olates the condition (27), rising to cTerr — 0.4 for the last 
epoch. 

3.4 Simulation error and Layzer- Irvine energy 

A test for the global accuracy of the integration of the equa- 
tions of motion is to measure how well the Layzer-Irvine 
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energy equation (Peebles 1980, equation 24.7) is obeyed (Ef- 
stathiou et al. 1985). One way to characterize this is through 
the change in the Layzer-Irvine integral, /, divided by the 
total potential energy W (Couchman et al 1995): 



I = K + W + j [2K + W\^ 



(28) 



where K is the total kinetic energy. In Tables 1 and 2 we 
present the percentage error in each of the simulations. The 
accuracy of the integration decreases as the spectral index 
steepens and as Q. decreases, the least accurate integration 
being that of the open n = simulation, for which the global 
error at the final epoch was of the order 1.8%. 



4 MEASURING THE POWER SPECTRUM 

The Fourier modes of the particle distribution can be deter- 
mined exactly using the expression (Peebles 1980) 



5k 



(29) 



Owing to the periodic boundary conditions, wavenumbers 
are restricted to be integer multiples of the fundamental 

mode, with an upper limit imposed by the finite sampling 
of the mesh: the Nyquist frequency. 



k 



Ny 



TT 

Ax 



(30) 



where Aa; = L/Njn is the mesh spacing and iVm is the dimen- 
sion of the mesh. The power spectrum can then be estimated 
through averaging over all of the modes in a thin shell in k 
space: 



m 

p{k)^{\s,f) = ^j2\^,,\ 



(31) 



where m is the number of modes to be averaged. This 
method is computationally inefficient, with the required cpu 
time scaling as MN for M modes. A faster method is to 
distribute particles onto a cubical mesh and perform a Fast 
Fourier Transform (FFT) . However, the assignment of mass 
to grid cells introduces some systematic effects which must 
bo corrected; these issues will be discussed in detail in Sec- 
tion 4.2. 

In the case of a 3D particle distribution, the task soon 
becomes memory limited, a 512"* FFT requiring roughly 
0.5 Gbyte for a 'real-real' transform and 1 Gbyte for a 
'complex-complex' transform. A solution to this problem 
was proposed by J98; we now detail this method, since it 
is critical for the present paper. 



4.1 Chaining the power 

Consider a ID discrete density field 6{x), which is periodic 

over a length scale L and which has a discrete Fourier trans- 
form given by equation (29) . If wo partition the density field 
using a coarse mesh with M grid colls, then the density at 
the point x can be described by the relation 



where x' is the position of the particle in its grid cell and j 
labels the cell. If we now map all of the grid cells into one 
cell, then the reduced density field, which is now periodic on 
the scale L/M, is 



5'-{x') = ^ 5{x' + jL/M) . 



(33) 



j=0 



The discrete Fourier transform of this reduced density field 
is then, 

N N 

^k = j^Yl exp(ifea;;) = j^Yl ""MiH^i - jL/M)] . (34) 

Provided that the fe-modes are integer multiples of the new 
fundamental mode, k = i2-w/{L/M), then the last term in 
the exponential is a multiple of 2n, so the modes of the 
reduced field are equivalent to the modes of the true field. 
There is, however, a reduction in the number of available 
modes, since the smaller volume of the coarse mesh gives a 
lower density of states. 



4.2 Numerical effects on the power 

There arc three important numerical effects which can mod- 
ify the 'observed' power spectrum from the true nonlinear 
signal: discreteness effects, charge assignment and force soft- 
ening. 



4-2.1 Discreteness effects 

For a random distribution of particles with no imposed clus- 
tering, the power does not vanish. This result can be deduced 
by splitting 3D space into a large rmmber of cubical cells, so 
that the occupation number of each cell is either rii = or 1 
(Peebles 1980). On computing the expectation of the power 
spectrum, we obtain the shot-noise spectrum 



which in dimcnsionless form is written 



- 

^«Vir.t — -rrr 



4tt 

TV 



kb 



(35) 



(36) 



6{x) = 6{x' + jL/M) 



(32) 



This leads us to write the true power spectrum, in the limit 
of large N (Peacock & Nicholson 1991) 

ALe(fc) = Albsik) - A^hot, (37) 

where A.'^^s is the observed power from equation (31). 

However, this correction is invalid for the glass and grid 
starts discussed in section 3. To determine the appropriate 
correction for these schemes wo directly computed the power 
spectrum of the initial conditions and then used these empir- 
ical spectra to construct a simple correction model. In Fig. 3 
(bottom) we show the raw power spectrum of the glass parti- 
cle load for the initial conditions and two subsequent epochs 
from the n = — 2 simulation. The glass power spectrum is 
characterized by a two-power-law spectrum: on intermedi- 
ate scales the spectrum is steep, roughly the n = 4 'minimal 
slope' (see section 28 of Peebles 1980) and at smaller scales 
this breaks to a shot noise spectrum. Furthermore, the bot- 
tom panel of Fig. 3 shows that the discreteness spectrum 
does not appear to evolve; we can therefore use the initial 
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1000 




1000 



1000 



Figure 3. (Top) The glass-discrctcncss corrected (squares) and 
uncorrected power spectrum (stars) of the glass n = — 2 simula- 
tion at an epoch a = 0.55. We show the linear fluctuation spec- 
trum (dashed line), which demonstrates that the box scale mode 
is still linear, the nonlinear spectrum according to the scaling 
formula of PD96 (thin solid line), a shot noise spectrum (dot- 
dashed line) and our two power-law discreteness model outlined 
in the text (thick solid line). (Bottom) Three epochs from the 
early stages of the same n = —2 simulation. From bottom to top 
epochs arc a = 0.025 (squares), 0.1 (circles) and 0.3 (stars). This 
demonstrates that the discreteness spectrum docs not evolve and 
also that the linear spectrum has been correctly established early 
on. Again, the lines are as in the top panel, with the thick solid 
line representing our fit to the discreteness spectrum. 



conditions to determine a discreteness correction that can 
be applied to correct the observed power at all subsequent 
epochs. This correction can be modelled as a transition be- 
tween shot noise on small scales and the almost minimal 
spectrum on intermediate scales: 



^glass 



(A^hot)"'''^ + (ALn) 



-1/a 



(38) 

with best fitting 



where a = 0.3 and A^^^ = (Ak/kh) 
values A = 0.0062 and = 6.8. 

For the grid, or '(jnict' start, the issue of a discreteness 
correction is fairly subtle, since there is initially no power 



added to the distribution by particle placement except on 
the scales of the Nyquist frequency of the mesh. However, 
as the simulation evolves under gravity, the sparseness of 
particles on small scales forms a power spectrum similar to 
a shot noise term on those scales. At late times this can 
be remedied by subtracting the Poisson spectrum from the 
raw power, since the large- and intermediate-scale modes in 
the evolved distribution are of higher amplitude than the 
shot-noise spectrum. At early times, when the true power is 
of relatively low amplitude, this approach is incorrect. Wo 
avoid the problem by excluding points whose amplitudes are 
below the Poisson spectrum on the equivalent scale. 

Fig. 4 (top) compares the uncorrected power spectra 
for the glass and grid starts measured from the n = —2 
simulations at two epochs a = 0.4 and a — 0.55. We observe 
that the nonlinear loci defined by the data for these two 
simulations are consistent and show no memory of the initial 
particle load. The only noticeable discrepancy between the 
two simulations is the difference in large-scale power; this 
arises because the simulations are independent realizations. 
Fig. 4 (bottom) contrasts the discreteness-corrected spectra; 
this shows that consistent final results axe obtained through 
simulating with grid or glass initial conditions. 



4-2.2 Mass assignment 

The assignment of mass onto the FFT mesh produces a fi- 
nite sampling error of the true density field. This problem 
was investigated for power spectra by Baugh & Efstathiou 
(1994), who proposed that equation (37) for the true field 
should be modified to 



o(fe) 



w(y) 



y = k/kn 



(39) 



where w{y) is the Fourier transform of the mass assignment 
window function, A^igj,(fc) is the appropriate discreteness 
correction and km = 2n/Ax is the wave number associated 
with the inter-mesh spacing Ax. However, we believe that 
there is a small flaw in their method. Any discreteness cor- 
rection should be made subsequent to the correction due to 
mass assigmnerit, since the discreteness correction accounts 
for the representation of a continuous field with a point-like 
distribution. We therefore implement the correction as 



w{y) 



(40) 



Several schemes exist for transferring mass onto the 
Fourier mesh. The simplest scheme is nearest grid point 
(NGP), which assigns all of the mass to the closest mesh 
point. More sophisticated methods such as cloud-in-cell 
(CIC) and triangular-shaped-cloud (TSC) attempt to smear 
the mass across a number of mesh points. We have adopted 
the TSC scheme to assign particles to the mesh. How- 
ever, the detailed correction is unimportant when using the 
chained-power method of J98. Results at high k can be ob- 
tained either by making substantial binning corrections to 
the main FFT mesh, or by moving to a sub-mesh of higher 
resolution. In practice, we make this transition before the 
corrections from binning become significant. Finally, Baugh 
& Efstathiou showed that, even after correcting for the win- 
dow function, the power is affected by aliasing close to the 
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Figure 4. (Top) Comparison of discreteness uncorrected power 
spectra measured from the quiet start (squares) and glass start 
(stars) n = —2 simulations at epochs a = 0.4 and 0.55. Line styles 
are as in Fig. 3. (Bottom) Comparison of discreteness corrected 
power spectra for the same outputs from the two simulations. 



Nyquist frequency. Again, when following the method de- 
scribed in Section 4.1, aliasing errors can be avoided by only 
using modes that are a safe distance from the Nyquist fre- 
quency of a given (sub)mesh (a factor of 2, in practice). 



4-2.3 Force softening 

The softening of the Newtonian force in the PP part of the 
A'^-body calculation (described in Section 3) induces an error 
in the integration of particle trajectories for close pairs. By 
considering the fractional error in the softened force from the 
true Newtonian force, we can impose some constraints on the 
small-scale cutoff, below which numerical effects dominate 
the clustering in our simulations. For our splinc-kcrncl force 
softening, we expect numerical effects to suppress the true 
power on scales of a few times the softening length. This 
corresponds to k/kh ~ 1700. 

The simplest way to discriminate between the true non- 
linear solution and numerical artefacts is to use the self- 
similar evolution of the scale-free simulations. Since the nu- 




Figure 5. Nonlinear power plotted against linear power (points) 
for the four scale-free simulations. For clarity, the data have been 
separated from each other by one order of magnitude in the y- 
direction, with the n = data untranslated. To determine the 

linear power given a nonlinear data point, the appropriate lin- 
ear scale is required. In the HKLM method, this is found using 
the transformation fcL = [1 + A^l('^'nl)]^^^''''snl- The solid line 
represents the fitting formula for the Einstcin-de Sitter models 
presented in Appendix B; the dashed line represents the PD96 
fitting formula; the dotted lines are the fits using the formula of 
JMW95. 



mcrical features are of fixed comoving length, the true den- 
sity field will scale under the transformations that were de- 
scribed in Section (2.4), whereas the numerical effects do 
not. We provide evidence for this in section 5. 



5 NUMERICAL RESULTS 
5.1 Similarity solution 

Fig. 5 shows the data for the four scale-free models in the 
HKLM form: nonlinear power on the nonlinear scale plotted 
as a function of the linear power on the linear scale. For 
clarity, the data have been separated from each other by 
one order of magnitude in the y-direction, with the n = 
data untranslated. In order to determine the linear scale and 
power that correspond to a given nonlinear data point, we 
use the nonlinear scaling relation (13). Explicitly, given a 
nonlinear data point fcjNL, Af^L) its linear counterpart is 

3+n 



feiL = (1 + ALl)"'/' 



kiL 
ko 



(41) 



where ko is a time-dependent normalization wavenumber de- 
fined by A^(fco) = 1 and we have assumed an initial power- 
law power spectrum for this example. 
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Figure 6. Nonlineaj" slope of the power spectrum versus spectral 
index. Points are measured from the scale-free simulation power 
spectra. The solid line represents the stable-clustering prediction. 
The dot-dash lines correspond to the predictions of the halo model 
(Ma & Pry 2000b), with an assumed e = 0.4 and 0.0 < /3o < 1.0. 
The dotted lines correspond to the halo model prediction with 
/3o = 0.25 and 0.4 < e < 1.0. 



When plotted in this form the scaling nature of these 

models is apparent. The power spectra measured from mul- 
tiple epochs of the simulations precisely overlay to define a 
single locus for each of the spectral models considered. We 
confirm the observation of JMW95 and PD96 that differ- 
ent spectral models produce different amounts of nonlinear 
growth and that the more negative the spectral index the 
steeper the locus in this plane. Fig. 5 also shows that the 
n = — 2 simulations have produced a single, tightly defined 
locus. This was not observed in previous studies (see Fig. 1 of 
Jain 1997, Fig. 1 of PD96 and Fig. 7 of Jain & Bertschinger 
1998). This failure of scaling in earlier n = —2 results was 
probably attributable to saturation of the box-scale mode. 

The evolution in the data can be roughly broken down 
into three regimes, the linear, the quasi-linear and the non- 
linear. General observations made about these regimes are: 

(1) Linear: A^l < 1: the 'nonlinear' power for all of the 
models converges to the linear power. 

(2) Quasi- linear: 2 < Aq < 30: the slope of the /nl curves 
are steep. Modelling the data in this regime with a single 
power-law of the form, Aq oc [Al]", we find for n = —2, 
— 1.5, —1 and 0, that the spectral slopes are a = 3.62±0.03, 
3.38±0.05, 3.12±0.06 and 2.96±0.1. This is reasonably close 
to the suggestion of Padmanabhan (1996) that Aq oc [Al]^, 
although there is a clear trend with n that is not expected 
in Padmanabhan's argument. This departure from a simple 
scaling relation is also supported by the results from loop- 
correction perturbation theory (see Fig. 19 of Scocciniarro 
& Frieman 1996). However, it may be argued that extended 
perturbation theory will fail at such large nonlinearities. One 
caveat is that it has been suggested that the nonlinear scal- 
ing relation may only truly be valid for ^ (Kanekar & Pad- 
manabhan 2001) and not A^. The small scatter observed in 
Fig. 5 leads us to believe that this might not be the case. 



(3) Nonlinear: A^l > 30: the /nl curves break away 
from the steep evolution which characterized the quasi-linear 
growth to form loci that are much shallower. Again, we have 
performed a simple power-law fit to the data of each locus. 

We find that for A'i,-, > 50 the n = — 2 data have a nonlinear 



slope a = 1.05 ± 0.09, and for Ai^ > 100 the n - 



-1.5, -1 



and data have nonlinear asymptotes of a = 0.87 ± 0.04, 
a = 1.08 ± 0.04 and a = 0.99 ± 0.04. This result is interest- 
ing for two reasons. Firstly, within the scatter in the simu- 
lations there appears to be little dependence on the initial 
spectrum for the nonlinear slope. Secondly, it is in clear con- 
tradiction to stable clustering, which predicts that a — 3/2. 
We note that this result agrees with the findings of Bagla, 
Engineer & Padmanabhan (1998) for clustering in 2D. How- 
ever, Fukusliige & Suto (2001) found that the stability on 
small scales, as measured from peculiar velocities, was not 
preserved locally but did apply globally. Our results do not 
agree with this. 

The shallow slope at high k may be interpreted in terms 
of the halo model. Ma & Fry (2000b) derived the following 
asymptotic limit for the power spectrum: 

18/3 - e(n 4- 3) 



7 = 



(42) 



2(3/3 + 1) ' 

where (3 ~ 0.8/3o is the power-law that governs the mass 

dependence of halo concentrations: c = r\./rs = {M / M,)'^" ; 
Vv and Ta being the virial and characteristic radius; and e 
is the power-law index that governs the low-mass tail of the 
mass function: dn/dM o<i v^;v oc l/a{M). Realistic values 
for e and /3o are 0.4 < a < 1.0 and 0.0 < /3o < 0.5. This is 
illustrated in Fig. 6, which shows the nonlinear power spec- 
tral index 7 as a function of the initial spectral index n. The 
values of 7 were obtained from the above nonlinear scaling 
relations, A^l oc [Al]", using the relationship (PD96) 

2>a{n + 3) 



Anl oc /Cnl! 



7 = 



(43) 



3 -I- a(n -1-3) 

We find 7 = 0.77, 0.91, 1.26, 1.49 for spectral indices 
n = — 2, —1.5, —1.0, 0.0. Comparing these measured values 
against the two predictions from equation (7) and (42), we 
see that 7 increases with the steepness of the spectrum, but 
that the data fall below the stable clustering prediction. In 
terms of the halo model, if one assumes e = 0.4 in accord 
with Sheth & Tormen (1999), then a strong dependence of 
/3o on n is required in order to match the measured data. 
On the other hand, if one adopts a value f3a = 0.25 in the 
middle of the current measured values, then it is impossi- 
ble to match the measured data with any value of e in the 
plausible range 0.4 to 1.0. In summary, equation (42) seems 
unable to predict the observed trend of 7(n) in a natural 
manner. This is puzzling, since we will show below that the 
general ideas of the halo model work very well in describing 
our data. One possibility is that equation (42) is valid only 
on scales smaller than those probed by current simulations. 

Also, in Fig. 5 we contrast our data with the fitting 
formula of JMW95 and PD96 (see Appendix Al and A2 
for these formulae). Both models work reasonably well in 
the quasi-linear regime, but with significant discrepancies. 
The n — —2 results are poorly fit by both models, with 
the power being in general underestimated; PD96 gives the 
poorer fit, and underestimates the power by up to a factor 2. 
The n = —1.5 locus is fairly well characterized by the JMW 
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function, but underestimated by PD96. The n — —1 results 
are fairly well fit by both models, except around the break 
between linear and quasi-linear slopes, where the functions 
overestimate the power. Finally, the n = locus is slightly 
overestimated at the linear to quasi-linear break by PD96 
and underestimated by JMW95. Wc have produced a new 
HKLM fitting formula that accurately fits the individual 
Einstein-de Sitter models, the results of which are shown 
in Fig. 5 as the thin solid line. The formula is described in 
Appendix B. 

5.2 Low-density power-law models 

In Fig. 7 we show how the nonlinear behaviour of the power- 
law models deviates from the scale- free solutions (solid lines) 
as the background density is lowered. Again, for clarity, the 
data have been separated from each other by one order of 
magnitude in the jy-diroction, with the n = data untrans- 
lated. In the linear regime, we again find that the nonlin- 
ear data follow the linear power. In the quasi-lineax regime, 
2 < Anl < 80, as O decreases, the locus defined by the 
data increases in amplitude relative to the scale-free models 
and the power-law slope steepens. This density-dependent 
evolution of /nl in the quEisi-linear regime was not appar- 
ent in previous studies (see PD96). The quasi-linear slope 
steepens as both n and f2 decrease. In the nonlinear regime, 
> 80, wc again observe that the slope of /nl is lower 
than the 3/2 value that is required by stable clustering. 

In Fig. 7, we also compare the data with the density 
dependent fitting formula of PD96. Again, the formula un- 
derestimates the shallow spectra and slightly overestimates 
the steeper spectra. However, the more striking discrepancy 
is that the formula suppresses the onset of density depen- 
dent growth until evolution is far into the nonlinear regime, 
and then tends to overestimate the highly nonlinear power. 
These discrepancies can in fact be seen in the comparison 
with the simulation data used by PD96. However, this li- 
brary of small (TV = 80^) simulations was in most cases 
unable to probe beyond A^l — 200, and so the deviations 
never became substantial. 

The failure of the JMW95 and PD96 functions to ac- 
curately model the Einstein-de Sitter data and account for 
the density dependence of nonlinear growth has clearly been 
shown. On attempting to fit this data set using the stan- 
dard HKLM-PD96 procedure we were able to produce an 
improved formula with an rms precision of 12%. However, 
on attempting to integrate the CDM models into the for- 
mulation, we could not find a satisfactory way to assign an 
effective spectral index to the models. We therefore decided 
to pursue an alternative approach to the problem of general 
nonlinear fitting functions, which proved to be more accu- 
rate. 



6 THE HALO MODEL FITTING FUNCTION 

In this Section, we attempt to describe the above nonlin- 
ear results by means of concepts abstracted from the 'halo 
model' (Peacock & Smith 2000; Scljak 2000; Ma & Fry 
2000a). The basic approach suggested by the halo model is 
to decompose the density field into a distribution of isolated 
haloes. Correlations in the field then arise on large scales 




Figure 7. (Top) HKLM plot for the open models. (Bottom) 
HKLM plot for the flat low-density models with a cosmological 
constant. Again, for clarity, the data have been separated from 
each other by one order of magnitude in the y-direction, with 

the n = data untranslated. For each model five epochs arc 
shown, these are: a = 0.6,0.7,0.8,0.9, 1.0, with the lowest locus 
for each model corresponding to the a = 0.6 epoch. In terms of the 
mass density parameter of tfic universe, these epochs correspond 
to: n = 0.294,0.263,0.238,0.217,0.200 for the open models and 
n = 0.619, 0.505, 0.407, 0.325, 0.260 for the A models. As in Fig. 
5 the solid lines represent the fitting formula for the Einstein-de 
Sitter models presented in Appendix B; the dashed lines represent 
fits from the PD96 function. 
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Figure 8. Comparison of tlic full lialo-model calculation as de- 
scribed in the text (thick solid lines) with the CDM data (points) . 
Also shown is the halo-model fitting formula from this work (thin 
solid lines). For clarity the four CDM models have been separated 
from each other by one order of magnitude in the y-direction, with 
the rCDM data untranslated. 



through the clustering of haloes with respect to each other 
and on small scales through the clustering of dark matter 
particles within the same halo. This then leads to a total 

nonlinear power spectrum 



PNL(fc) =PQ(fc)+PH(fc), 



(44) 



where PQ{k) is the quasi-linear term that represents the 
power generated by the large-scale placement of haloes and 
where Pnik) describes the power that results from the self- 
correlation of haloes. 

It is remarkable that such a simple decomposition ap- 
pears to work well in describing the main characteristics of 
the two-point correlations of the cosmological mass density. 
It is maybe yet still more impressive when one considers that 
the present formulation knows nothing of the large-scale fil- 
amentary structure of the density field (which is governed 
by the correlation function of halo centres) . Indeed this de- 
ficiency was recently pointed out and addressed by Scocci- 
marro & Sheth (2002) 

In Fig. 8 we directly compare the halo-model calcula- 
tions (thick solid lines) with the CDM simulations of J98 
(data points - see section 6.4 for full description). Also 
shown is the halo-model fitting function that we present 
later (thin solid lines - see Appendix C). The halo model 
calculations are exactly those of Peacock & Smith (2000). 
Prom the figure, it can clearly be seen that the calculations 
qualitatively reproduce the data for all of the models, but 



that in detail only match SCDM and rCDM closely. Fur- 
thermore when one attempts to model the power-law spectra 
the results are worse, with the n = case being an extreme 
example (see the later discussion). Thus our aim in what 
follows will therefore be to produce a simple fitting formula 
that draws on the broad elements of the halo model, such 
as the above decomposition of the power spectrum into two 
linearly summed terms, but which is of very high accuracy. 



6.1 The quasi-linear term 

Consider the quasi-linear part first. Seljak (2000), Ma & 
Fry (2000a) and Scoccimarro ct al. (2001) assumed that 
one should use linear theory filtered by the effective window 
corresponding to the distribution of halo masses, convolved 
with their density profiles and a prescription for their bias 
with respect to the underlying mass field: 



dM h{M) n{M) p{k, M) 



(45) 



where n (M) dM is the mass function, p (fc, M) is the Fourier 
transform of the density profile and b{M) is the bias field of 
dark matter halo seeds. Peacock & Smith (2000) made the 
simpler assumption that the quasilinear term corresponded 
to pure linear theory: 



PqW = PL(fc). 



(46) 



This is equivalent to equation (45) on large scales, since in 
this limit the filtering effect of haloes is negligible, and we 

must have 



dM b{M)Mn{M) = 1 



(47) 



Neither of these approaches is really satisfactory, since 
Ph comes to dominate only at scales where linear theory 
must break down to some extent (Al ~ 1). Quasi-linear 
effects must modify the relative correlations of haloes away 
from linear theory, irrespective of whatever allowance may 
be made for the finite sizes of haloes. One way of seeing this 
is via the scaling part of the HKLM procedure: see equations 
(13). This shift of scales from gravitational collapse causes 
a significant change in power at wavenumbers where Al is 
of order unity - which is just the point where the filtering 
effects of the largest haloes will also start to bo important. 
An alternative point of view is provided by perturbation 
theory, which suggests that quasilinear effects should tend 
to suppress power for n > —1.4, but enhance power for 
more negative indices (e.g. section 4.2.2 of Bernardeau et 
al. 2001). Again, such effects cannot be cleanly separated 
from the convolving effects of halo profiles. We therefore 
take an empirical approach, allowing the quasilinear effects 
to depend on n. Since the philosophy of the halo model is 
that Aq should be negligible on small scales, we also build 
in a truncation at high k: 



l + a„A2(fc) 



exp-f{y)- y. 



(48) 



where k^ is a nonlinear wavenumber, defined below in Sec- 
tion 6.3 a„ and /3„ are spectral dependent coefficients and 
f{y) is the polynomial y/4 + y'^/8 that governs the decay 
rate. We adopt this expression for all spectra. 
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6.2 The halo term 

In the halo model the self-halo term is (Seljak 2000; Peacock 
& Smith 2000; Ma & Pry 2000a; Scoccimarro et al. 2001) 



p2 (2^)^ 



dM n(M) \p{k,M)f 



(49) 



In order to model this we want something that looks like 
a shot-noise spectrum on large scales, but is progressively 
reduced on small scales by the filtering effects of halo profiles 
and the mass function. In terms of the dimensionless power 
spectrum, a candidate form for this is 



an y 



1 + 6„ y + C„ t/3-7„ 



y = k/ka 



(50) 



jn ) are dimensionless numbers that de- 
pend on the spectrum. However, with Pn defined in this 
way, the formalism defined by equation (44) breaks down 
for steep spectra. The self-halo power clearly dominates at 
small k for any spectrum that is asymptotically r?, > (e.g. 
all CDM models). This has been independently noted by 
Sheth & Cooray (2001). The halo model thus fails to re- 
spect low-order perturbation theory in such cases, and this 
is a clear defect of the model. 

In order to solve this problem, the self-halo power must 
become steeper than Poisson on the largest scales. This 
malces sense if we think of the halo model as a two-stage 
process: (i) fragment a uniform mass distribution into a set 
of haloes; (ii) move these haloes according to a superimposed 
large-scale displacement field. Since the first stage conserves 
mass, the large-scale power spectrum must approach a 'min- 
imal' form with n = 2 (e.g. section 28 of Peebles 1980). If one 
conserves momentum also, the minimal spectrum becomes 
even steeper: n = 4. It is a moot point which of these is the 
appropriate asymptote for this problem, since the two-stage 
view of the halo model is only a heuristic argument. Since 
we will never wish to consider spectra that are asymptoti- 
cally much steeper than n = 1, it will suffice to force the 
n = self-halo term to approach n = 2 on sufficiently large 
scales. This can be achieved if equation (50) is modified as 
follows 



A|(fe) 



^l'{k) 



1 + /Un2/"^ +Vny~ 



y = k/ka 



(51) 



where we have introduced a term in fc* in order to soften the 
transition to the k^ slope. Again, the parameters /u„ and Vn 
are spectral dependent coefficients. 



6.3 The nonlinear scale 

In order to implement these arguments, we need an appro- 
priate general definition of the nonlinear scale (see Section 
2.4), which should be related to the characteristic mass in 
the halo mass function. As studies over many years have 
shown with increasing accuracy (Press & Schechter 1974; 
Sheth & Tormcn 1999; .Icnkins et al. 2001), the halo mass 
function appears to depend only on the dimensionless fluc- 
tuation amplitude 



(52) 



where 5c is a constant of order unity, usually identified with 
the linear over-density for collapse in the spherical model 



Table 3. The cosmological parameters of the N = 256^ CDM 
simulations from J98. For these CDM models F = Q.h is the shape 
parameter of the spectrum, its is the rms fluctuation in spheres 
of 8/i~^ Mpc and h is the Hubble parameter 



Model 


F 


0-8 


il 


A 


h 


L//i-iMpc 


SCDM 


0.50 


0.51 


1.0 


0.0 


0.5 


239.5 


SCDM 


0.50 


0.6 


1.0 


0.0 


0.5 


84.55 


rCDM 


0.21 


0.51 


1.0 


0.0 


0.5 


239.5 


tCDM 


0.21 


0.6 


1.0 


0.0 


0.5 


84.55 


ACDM 


0.21 


0.90 


0.3 


0.7 


0.7 


239.5 


ACDM 


0.21 


0.90 


0.3 


0.7 


0.7 


141.3 


OCDM 


0.21 


0.85 


0.3 


0.0 


0.7 


239.5 


OCDM 


0.21 


0.85 


0.3 


0.0 


0.7 


141.3 



and R is the effective filter radius. The multiplicity function 
for haloes, which is defined as the fraction of mass carried by 

haloes with mass in a logarithmic interval, peaks for systems 
where a(R,t) is of order unity, and we can therefore choose 
to define the nonlinear scale in this way: 



a{k-\t) = l. 



(53) 



This definition of scale depends on the functional form cho- 
sen to filter the spectrum, but the main effects of changes in 
this choice can be absorbed into the fitting coefficients. We 
therefore take the convenient choice of a Gaussian filter: 



/ 



cT^(i?G,t)= / Ai(fc,t) exp(-A;^i?|,) dlnfc. 



(54) 



With this choice of filter, scale-free spectra have 

3+n 

A^(fc,t) 

-l/(3-Hn) 



k 

ko{t) ) 

[{l+n)/2\\ 



(55) 



6.4 Application to CDM 

We have generalized our formula to fit the Virgo and GIF 
CDM simulations from J98, which comprise four models: 
SCDM; tCDM; ACDM; OCDM. Table 3 lists the cosmo- 
logical parameters for these models. The data are publicly 
available from http://www.mpa-garching.mpg.de/Virgo/. 
We have re-measured the power spectrum for the epochs 
z = 0.0, 0.5 1.0 2.0 and 3.0 for both the Virgo and GIF data, 
the results are presented in Figs. 14 and 15. The transfer 
function for these simulations was that of Efstathiou, Bond 
& White (1992): 



A^(fe) = 



Ak" 



{l + [ag+(6g)3/2 + (cg)2].}2/. 



(56) 



where q = k/V, a = 6.4 /i"^Mpc, b = 3/i"^Mpc, c = 
1.7 Mpc. The normalization constant A is chosen by fix- 
ing 0-8 . 

In order to model these more general curved spectra, 
we define an effective spectral index via 



3 -I- neff = - 



din a'^{R,t) 



din J? 



(57) 
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Figure 9. Variation of effective spectral index (top panel) and 
curvature (bottom panel) as a function of the rms fluctuation in 
Gaussian spheres of radius Rq, for the four cosmological models 
considered. Note the lower erg values that corresponded to the 
big-box simulations has been assumed for the SCDM and tCDM 
models. 



Since the mass function should depend mainly on the Tay- 
lor expansion of o about the nonUnear scale, we also allow 
dependence on the spectral curvature: 



C ; 



(58) 



For the case of a Gaussian filter these expressions have the 
explicit forms, 



3-|-neff = ^ / dlnfc AL(fc,t) y^exp (-j/^'' 



and 

C = {Z + n,sf 



(59) 




r=nh 

Effective spectral index n^jj 



I 




I 



a{Ry 



Figure 10. Dependence of the nonlinear wave-number ka (top 
panel), effective spectral index n^f[ (middle panel) and curvature 
parameter C (bottom panel) on the shape parameter F and nor- 
malization erg of the linear power spectrum. The parameters n^g 
and C are degenerate under F and erg. This degeneracy is, how- 
ever, broken by the nonlinear wavenumber. 
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Table 4. The nonlinear wavenumber ka in units of h Mpc"'^, the 
effective spectral index rieff and curvature C of the spectrum on 
the nonlinear scale, for the four CDM models listed in the text. 



Model 






C 


SCDM 


0.574 


-1.455 


0.411 


rCDM 


0.735 


-1.850 


0.305 


ACDM 


0.306 


-1.550 


0.384 


OCDM 


0.332 


-1.581 


0.375 



j d\nk Ai{k,t) (j/^ - j/"*) exp (-j/^) 



(60) 



where y — kRc and where the exphcit time dependence of 
the power spectrum has been kept to indicate the redshift 
dependence of the elfective quantities. In Table 4 we list 
the nonlinear wavenumber, effective spectral index and cur- 
vature of the spectrum on the nonlinear scale for the four 
Virgo (big-box) CDM models, generated according to the 
above prescription. 

Fig. 9 shows the variation of the effective spectral index 
(top panel) and curvature (bottom panel) for the four Virgo 
CDM models with the rms fluctuation measured in Gaussian 
spheres of effective radius Rq- The effective spectral index 
is quite sensitive to whether it is defined at tr = 1 or at 
some other value. However, including the curvature (which 
depends much more weakly on a), means that this uncer- 
tainty is automatically allowed for. With the nonlinear scale 
and effective spectral index and curvature as defined through 
equations (53-58), we find that we can accurately model 
CDM spectra. As expected, Fig. 9 shows that the OCDM 
and ACDM models are almost indistinguishable: both pos- 
sess nearly identical linear power spectra, with only a slight 
difference in normalization. The tCDM model has the shal- 
lowest effective spectral index, almost approaching n = —2 
and the SCDM model has the steepest, with n = —1.4. 
The power-law models that we have simulated encompass 
this range of ries. Thus, we are confident that the new fit- 
ting function will be constrained by the appropriate range 
of spectral models, with the notable exception of the z > 3 
tCDM data for which rics < —2. These models are the sole 
basis for the fitting formulae in the n < —2 regime. 

Fig. 10 shows the dependence of ka (top panel), riefF 
(middle panel) and C (bottom panel) on the shape param- 
eter and normalization of the linear power spectrum for 
(0.1 < r < 0.8) and (0.4 < ag < 1.2). In aU of the models 
dark contrast represents a higher value. The parameters riefi 
and C are degenerate under T and as. This degeneracy is, 
however, broken by including the nonlinear wavenumber. 



6.5 Parameter optimization 

We now give the best-fitting coefficients, including depen- 
dence on cosmology. These coefficients were obtained by op- 
timizing the formula to fit the scale-free and Q < 1 power- 
law simulations described here; the CDM simulations of J98; 
and on large scales {k < 0.15ftMpc^^). the results of 2nd- 
order perturbation theory (calculated using the formulae of 
Lokas et al. 1996). Owing to the fact that numerical simula- 
tions are susceptible to sample variance on large scales, ana/- 
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Figure 11. (Top) nonlinear power ratioed to the linear power as 
a function of wavenumber scaled in terms of the normalization 
wavenumber Aiq: where A^(fco) = 1- The data points are for the 
scale-free simulations; the solid lines represent the fits from the 
new halo based formula in Section 6.5; the dotted lines are PD96 
fits. (Bottom) The goodness of the new fit. The y axis represents 
the ratio of observed nonlinear power to nonlinear power pre- 
dicted by the halo based fitting function. The x axis is observed 
nonlinear power. 



lytic perturbation theory results were preferred. In the halo 
model the cosmology dependence arises in a subtle way. To 
the extent that the mass function depends only on v (when 
expressed as a function of -R) and that 6c has no strong 
cosmology dependence, the mass function for a given spec- 
trum is also independent of cosmology. Therefore, the only 
effect on the halo power spectrum should be through the 
sizes of haloes; these depend on cosmology because haloes 
that collapse at high redshift are smaller. Collapse redshift 
is a function of mass and cosmology (see e.g. Appendix C 
of Peacock & Smith 2000). High-mass haloes always have 
Zc ~ 0; these thus filter the large-scale part of the spec- 
trum in a cosmology-independent way. Conversely, low-mass 
haloes are important at high k, and these do depend on cos- 
mology - which alters the effective scale at which filtering 
occurs. However, there appears to be no simple way to im- 
plement such a complicated dependence into the fitting pro- 
cedure. We therefore insert empirical functions of fl into the 
procedure. Also, motivated by the findings of Section 5, we 
allowed the power-law indices that govern the quasi-linear 
regime to be density dependent. 
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Figure 12. Top and bottom panels arc similar to Fig. 11, but tiiis 
time points represent open model data. Five epochs are shown; 
these arc a = 0.6, 0.7, 0.8, 0.9, 1.0. In terms of Q, these epochs cor- 
respond to: n = 0.294,0.263,0.238,0.217,0.200. The thick solid 
line represents the new halo-model based fitting to the scale-free 
data. 
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Figure 13. Top and bottom panels arc similar to Fig. 11, but 
this time points represent A model data. Four epochs are shown; 
these arc a = 0.7,0.8,0.9,1.0. In terms of Q, these epochs cor- 
respond to: n = 0.505, 0.407, 0.325, 0.260. Again, the thick sohd 
line represents the fit to the scale-free data. 



The prescription that was found to work best is given 
in Appendix C. Code to evaluate the fitting function can 
be downloaded from the web address listed in the abstract. 
Note that the above coefficients were obtained by fitting the 
data over a restricted range of scales. The scale-free data 
were constrained to have k/ka > 0.3. The open and A data 
were constrained to lie in the range : (4.0 < A?, < 15.0) for 
n = -2; (0.3 < A?, < 15.0) for n = -1.5; (0.3 < A?. < 20.0) 
for n = -1; (0.3 < A?, < 25.0) for n = 0. The CDM data 
were fit under the constraints: k > 0.3 for the big box data 
and to > 0.5/iMpc~'^ for the higher resolution small box 
calculations; the nonlinear power must be 10% greater than 
the discreteness correction, equation (38). On larger scales 
k < 0.15/iMpc~^, the formula was calibrated to the results 
of 2nd-order perturbation theory. 

In Fig. 11 we compare our new halo-based fitting func- 
tion with the scale-free simulations. The new model clearly 
reproduces the data to a high degree of accuracy. Also, it 
is important to note that when the data are plotted in this 
way the scaling nature is again apparent and the departure 
from stable clustering, which is indicated by the deviation 
away from PD96 for k/ko > 10, is pronounced. 

In Figs 12 and 13 we compare the new halo based model 
with the power-law data for O < 1 and fi -|- A = 1. For 



all of the models the inclusion of the functions /i, /2 and 

/a, seems to well reproduce the observed density-dependent 
growth. The only significant discrepancy is for the n = —2 
open data, where the power is underpredicted in the quasi- 
linear regime. 

In Figs 14 and 15 we compare the model with the CDM 
data. Again, the model does exceptionally well at repro- 
ducing all of the data over the range of scales where we are 
confident that numerical effects are unimportant. In particu- 
lar, the OCDM and rCDM predictions are very significantly 
improved using the new prescription. 

Having demonstrated the success of the halo fitting 
function on small scales, we next consider the large scales. 
We assess this using the predictions derived from 2nd order 
perturbation theory (see Baugh & Efstathiou 1994). Fig. 16 
shows the ratio of nonlinear to linear power for four CDM 
models. The current models match perturbation theory for 
k < O.l/iMpc"^, but deviations exist at higher k. These 
plausibly reflect a genuine breakdown of perturbation the- 
ory, since the model was required to match perturbation 
theory as well as possible for k < 0.15 ftMpc^^, and yet the 
fit is breaking down slightly before this upper limit. Both 
the halo fitting function and PD96 agree well in this range. 
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Figure 14. Power spectra for the four Virgo CDM simulations (J98) in large cosmological volumes, L = 239.5/i~^Mpc. Each panel 
shows the evolution of structure with redshift. The data points correspond, from low to high, to epochs 2 = 0, 0.5, 1.0, 2.0 and 3.0. Note 
that only those points with a measured power above the discreteness spectrum are plotted. The solid line represents the new halo-model 
based fitting procedure, with dotted lines representing the decomposition into the self-halo and halo-halo terms; the dashed line is the 
PD96 fit. 



7 CONCLUSIONS AND DISCUSSION 

In this paper we have presented a set of high-resolution, 256^ 
particle, scale-free A'^-body simulations, designed to investi- 
gate self-similar gravitational clustering and in particular 
the effects of nonhnoar evolution. Wc have also performed a 
further series of numerical simulations, with the same res- 
olution, to explore how the evolution of clustering depends 
upon the background density of the universe. Together, these 
simulations represent the best calculations that exist to date 
for the set of models explored, with a factor 512 improve- 
ment in mass resolution over the ground-breaking work of 
Efstathiou et al. (1988). 

We verified that the final output power spectra wore ro- 
bust by considering grid and glass particle loads. However, 
at early times the problem of discreteness correction is sim- 
pler to handle if a glass start is applied; we have described 



a detailed method for correcting the clustering signal in this 
case. We have implemented the power spectrum estimation 
technique of J98, which allowed us to probe high spatial 

frequencies without aliasing effects or errors due to mass as- 
signment to the Fourier mesh. The simulation results may 
be summarized as follows: 

(1) Scale-free simulations with < n < — 2 show self- 
similarity under the scaling fco(a) oc a"^/^""*"^'. This conclu- 
sion is in agreement with the results of Efstathiou (1988) 
and Jain & Bcrtschingor (1998). 

(2) In the quasi-linear regime, the power spectrum is 
characterized by a steep power law. The exact slope depends 
upon the spectral index n of the input spcctrxim and the 
value of O, the slope steepening as n becomes more negative 
and as O is reduced. 

(3) The observed nonlinear asymptote of the Einstein- 
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Figure 15. Same as for Fig. 14, but tliis time for the smaller box GIF simulations. 



de Sitter simulations was found to bo inconsistent with the 
Anl [Al]^^^ prediction of stable clustering. A shallower 
slope with Anl oc [Al]^ is preferred. This result makes 
sense in terms of the halo model: calculations using the ex- 
tended Press-Schechter apparatus show that haloes will tend 
to merge with systems of similar mass to their own (Laccy & 
Cole 1993). Mergers of this kind will disrupt the virial equi- 
librium of the system, violating the basic assumption that 
underlies stable clustering. However, if this process were rare 
then stable clustering could be upheld in a statistical sense. 

(4) The nonlinear fitting formulae of PD96 and JMW95 
failed to reproduce the n = —2 results and were only 
marginally successful at reproducing the steeper spectra. 
The low-density power-law data were poorly fit by PD96. 

(5) For the $1 < 1 simulations, it is interesting to con- 
sider how the nonlinear slope changes with density. In the 
nonlinear limit equation (C4) (appendix C) becomes 

A^{k) oc k . (61) 

For a given n, fi increases as Q decreases, and so the power- 



law slope steepens. This result supports the idea that small 
scale clustering is more closely related to the emergence of 
the internal density structure through the continual accre- 
tion and merger of haloes. The reasoning is as follows: for a 
low-density universe mergers are less frequent and so haloes 
have more time to virializc. This means that stable cluster- 
ing may be considered to be a better approximation for these 
systems. Prom the arguments in Section 1 and 5, this would 
then be manifest as a steepening of the nonlinear slope. 

In the second part of this paper, we proposed an im- 
proved fitting function for mass power spectra to replace 
the much-used PD96 formula. We have adopted a new ap- 
proach to fitting power spectra, based upon a fusion of the 
halo model and a HKLM scaling. The method was general- 
ized to fit more realistic curved spectra, by introducing two 
new parameters, rics the effective spectral index on the non- 
linear scale, and the spectral curvature, C. We found that 
the halo model as previously envisaged in the literature fails 
to approach linear theory on large scales for n > 0. We 
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Figure 16. Comparison to perturbation theory. The top panel 
shows the ratio of the evolved power to the linear power for two 
CDM models with F = 0.2, but with different normalizations. 
The model with as = 1. has been translated by a factor of 0.5 
in the j/-direction. The points represent perturbation theory; the 
thick solid line is this work; the dashed line represents PD96. The 
bottom panel is the same as the top, but with F = 0.5. 

have argued that this should be cured by changing the self- 
halo power from n = to n = 2 on large enough scales, 
and we have shown empirically that this approach allows 
an accurate description of a very wide range of power spec- 
trum data. Our new fitting formula reproduced the scale- 
free power spectrum data and also the CDM results of J98 
with an rms error better than 7%. This is to be preferred to 
the widely-used PD96 prescription, and should be useful for 
a variety of cosmological investigations. In particular, our 
preliminary investigations show that the present formalism 
should cope naturally with spectra containing a realistic de- 
gree of baryonic features (e.g. Meiksin, White & Peacock 
1999). 

The halo model provides a novel way to view structure 
formation, and has yielded useful insights into the origin of 



nonlinear aspects of galaxy clustering. This work heis con- 
centrated on the low-order statistics of the density field, but 
it is also possible to consider higher-order statistics such as 
the bispectrum. This three-point function in Fourier space 
probes the shapes of largo-scale structures that arc gener- 
ated by gravitational clustering. No shape information is 
included in the current formalism, so it will be interest- 
ing to see how well the model can account for higher-order 
statistics. Initial results in this direction (Scoccimarro et al. 
2001) seem to be promising. In general, the important ques- 
tion is the extent to which the halo model can encapsu- 
late the phase information in the density field, since fields 
with identical power spectra can possess completely differ- 
ent real-space distributions (e.g. Chiang & Coles 2000). The 
halo model will inevitably fail to encompass these details of 
the density field in full, although it may still offer useful in- 
sights. However, at the two-point level, we have shown that 
the model is far more than an educational device, and it 
can be used as a tool for a high-precision description of the 
evolution of the dark-matter power spectrum. 
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APPENDIX A: HKLM FITTING FUNCTIONS 
Al The JMW95 function 

The JMW95 function was designed to model the n- 

dcpcndcnce of the nonlinear evolution of scale-free power 
spectra. The formula was also used to model fl — 1, CDM- 
like models through the adoption of an effective spectral 
index; see equation (14). JMW's formula described their nu- 
merical data with an rms accuracy of 15-20%, but for our 
higher resolution scalc-froo data the fit is much worse, hav- 
ing an rms error of 56%. Their formula is 



:.(^nl) 



B{n) 



Ai(fcr 



B{n) 



(Al) 



where B{n) is a constant which depends upon the spectral 
index n and where /jmw (y) remains independent of n. The 
explicit forms are: 

3 + n\^-'^ 



B(n) 
and 



/jMw(y) = y 



where y = Al{kj^)/B{n) 



1 + 0.6y + y' - 0.2y' 



1 -I- 0.0037J/3 



(A2) 



(A3) 



A2 The PD96 function 

PD96 performed a similar study to JMW95, but extended 
the set of cosmological models to include f2 < 1 open and 
flat universes. They also improved on JMW95 by including 
CDM data in the optimization procedure and by proposing 
that the effective spectral index would vary continuously 
with scale: equation (15). They reported that their fitting 
formula described their simulation data to an accuracy of 
about 14%, but it describes our complete data set with an 
rms error of 54%. The PD96 fitting formula is 



/pd(j/) = y 



1 + Bl3y+ [Ay_ 



l + ([^j/]"g3(n,A)/[l/yi/2])/3 



1//3 



(A4) 



where y = Ai^{k^). B describes a second order deviation 
from linear growth; A and a parameterize the power law 
that dominates the function in the quasi-linear regime; V 
is the virialization parameter that gives the amplitude of 
the /nl oc y^^^ asymptote; /3 softens the transition between 
these regions; g{ii) is the density dependent growth factor 
of (Carroll, Press & Turner 1992), which is the ratio of the 
linear growth factor to the expansion factor. This has the 
functional form 

5" 1 
= ^n[n*/^-A-f (l + n/2)(l-FA/70)] . 



S(«) 



The best-fitting parameters were 

A = 0.482(1 -I- n/3)-°-^*^ 

B = 0.226(1 -I- n/3)"^ '^''* 

a = 3.310 (l-Fn/3)-°-^'''' 

= 0.862(1 4- n/3)"°-^*^ 

V = 11.55 (l-^r^/3)-°■^^^ 



(A6) 
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APPENDIX B: NEW HKLM FITS TO THE 
PRESENT DATA 

We have performed a nonlinear least squares fitting to the 
individual scale-free loci (see Fig. 5) using a single formula. 
The individual fitting functions are accurate to ~ 9%. The 
formula is 



/Eds(y) = y 



l+y/a + {y/bf + {y/cr 



1 + 



1/7 



(Bl) 



where y = AlC^Cl) and the relevant parameters for each n 
are presented below 



n 

-2 
-1.5 
-1 




a 

3.138 
2.710 
10.37 
29.26 



h 

0.358 
0.710 
1.115 
1.394 



c 

0.527 
0.919 
1.403 
1.941 



d 

0.940 
1.852 
2.873 
3.753 



a 
8.247 
0.707 
6.655 
6.547 



P 7 

0.508 0.330 

0.647 0.332 

0.697 0.366 

0.847 0.351 



APPENDIX C: THE HALO MODEL FITTING 
FUNCTION 

The halo model decomposes the power into a sum of two 

contributions: 



AL(fc) = A?j(fc) + A|(fc). 
These are given separately by 

{1 + Al{k)f- 



A%{k) = Al{k) 



exp[-/(j/)]; 



l + a„Al{k) 
where y = k/ka and f{y) = y/A + y^/8; and 
Al'(fc) 



A%{k) 



where 



1 + Mn2/"^ + i^ny~ 



an y 



3/i(0) 



l + 6„j//2(o) + [c„/3(n) y] 



3-7n 



(CI) 



(C2) 



(C3) 



(C4) 



and y = k/kc,. 

The parameters of the spectrum are defined via Gaus- 
sian filtering: 



Al{k) exp{-k'^Rl) dink. 



In these terms, 
a{k-') = 1. 
The effective index is 
din 0-2(7?) 



3 + ries 



dlnR 

and the spectral curvature is 



C 



din i?2 



(C5) 



(C6) 



(C7) 



(C8) 



Allowing (a„, bn,Cn,'yn,Oin, Pn, IJtn, v„) to Vary as a function 
of spectral properties, the following coefficients fit our sim- 
ulation data and the CDM simulations of J98 to an rms 
precision of 8.6% (very much better than PD96). In partic- 
ular, the model describes the ACDM data of J98 extremely 
well. For redshifts z < ?>, the deviation in power between 



model and the average of the large-box and small-box data 
from J98 is always less than 3% for k < lO/iMpc"^. This 
represents a perfect fit with present knowledge, since the two 
datasets themselves can differ by at least this much. Note 
the use of terms up to in the fit for a„; these are required 
in order to describe the rapid rise in amplitude of the halo 
term for n < —2. For less negative n, the higher-order terms 
arc unimportant. The coefficients are: 

logio «n = 1-4861 + 1.8369 n + 1.6762 + 0.7940 

+0.1670 n* - 0.6206 C ; (C9) 

logio = 0.9463 + 0.9466 n + 0.3084 - 0.9400 C ; (CIO) 

2 0.0793 C ;(C11) 

(C12) 



log^gCn = -0.2807-^0.6669 n-F0.3214 n 



7„ = 0.8649 + 0.2989 n + 0.1631 C : 
an = 1.3884 -I- 0.3700 n - 0.1452 
Pn = 0.8291 -I- 0.9854 n + 0.3401 
logio = -3.5442 -I- 0.1908 n ; 
log^o fn = 0.9589 + 1.2857 n ; 
and the dependent functions are: 

n < 1 



flam 


= n 


-0.0732 


/2a (f2) 


= Q 


-0.1423 


/3a (n) 


= n 


0.0725 




= n 


-0.0307 


f2b{n) 


= n 


-0.0585 


f3bm 


= n 


0.0743 



n-i-A = 1 



(C13) 
(C14) 
(C15) 
(C16) 

(C17) 
(C18) 



For models in which A is neither zero nor 1 — fl, we suggest 
interpolating the functions /i etc. linearly in A between the 
open and flat cases. 
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